Modeling the effects of increasing the testing capacity via a SIR-like model

In [1]:
import os
import git

# Our numerical workhorses
import numpy as np
import scipy as sp
import pandas as pd

# Import Interactive plot libraries
import bokeh.plotting
import bokeh.layouts
from bokeh.themes import Theme
import holoviews as hv
import hvplot
import hvplot.pandas
import panel as pn

# Import the project utils
import covid19

bokeh.io.output_notebook()
hv.extension('bokeh')
Loading BokehJS ...
In [2]:
# Set PBoC style for plot
theme = Theme(json=covid19.viz.pboc_style_bokeh())
bokeh.io.curdoc().theme = theme
hv.renderer('bokeh').theme = theme

Objective

In this notebook we'll explore a simple SIR-like model to better understand the effect that an increasing testing capacity would have in isolating the COVID-19 cases. Here we'll present interactive plots in which parameters can be modified via sliders. This interactive feature facilitates the exploration of the complicated parameter space.

The model

For this model we will use the well established SIR-like model in which the population is compartmentalized into the following groups:

  • $S$: Population susceptible to be infected by the SARS-CoV-2 virus.
  • $E$: Population that has been exposed to an infected person long enough to contract the disease. This can be thought as the incubation time between the contagious interaction and the infective phase.
  • $I_S^{(q)}$: Infected population that presents symptoms sever enough to be set in quarantine (subscript $(q)$ denotes quarantine).
  • $I_A$: Asymptomatic infected population that, unless tested, wouldn't know they are infected and spreading the virus.
  • $I_A^{(q)}$: Asymptomatic infected population that, after being tested positive for COVID-19, complies with the quarantine order.
  • $R$: Population that either recovered or unfortunately died due to the disease.

Note that we do not distinguish between people that recovered form the disease or that passed away. As far as the model goes, there is no distinction to how the dynamics of the epidemic change. Nevertheless if needed the model could easily divide the $R$ population into these two groups.

The assumption of our model are the same as the typical SIR model. What this means is that we assume we operate at the limit of large populations where small fluctuations are averaged out, so we use deterministic ODEs. We also assume a well-mixed population where contacts follow mass-action dynamics. All these assumptions can be easily relaxed, but for now we just want to get the general possible scenarios that can result from running the dynamics with different parameter regimes.

Let's now write down the equations for each of the populations. First for the susceptible people we have

$$ {dS \over dt} = - \beta_S S I_S^{(q)} - \beta_A S I_A - \beta_A^{(q)} I_A^{(q)}, $$

where $\beta_s$, $\beta_A$, and $\beta_A^{(q)}$ are the infective rates for each of the infective populations.

$$ {dE \over dt} = \beta_S S I_S^{(q)} + \beta_A S I_A + \beta_A^{(q)} I_A^{(q)} - \gamma_S E - \gamma_A E, $$

where $\gamma_S$ and $\gamma_A$ are the rates at which exposed people become infective presenting symptoms or being asymptomatic, respectively. For the three infected populations we have

$$ {dI_S^{(q)} \over dt} = \gamma_S E - \mu I_S, $$$$ {dI_A \over dt} = \gamma_S E - \theta I_A - \mu I_S, $$$$ {dI_A^{(q)} \over dt} = \theta I_A - \mu I_S, $$

where $\theta$ is the rate in which an asymptomatic person goes into quarantine. This rate will depend on the availability of tests, the accuracy of such tests, and the probability that a person actually complies with the quarantine. We also assume that all infected groups either recover or die at the same rate $\mu$. With this the dynamics for the removed group are of the form

$$ {dR \over dt} = \mu \left(I_S^{(q)} + I_A + I_A^{(q)} \right) $$

To look at the dynamics we will numerically integrate the system of ODEs using scipy's ode_int function. For this we first need to define a function that computes the right-hand side (rhs) of

In [3]:
def system_rhs(x, t, beta_s, beta_a, beta_aq, gamma_s, gamma_a, theta, mu):
    """
    Right hand side of system. Returns time derivatives
    Parameters
    ----------
    x : array-like. len(x) = 6
        Array containing all 6 compartmentalized populations in the following 
        order: S, E, I_s, I_a, I_aq, R.
    t : array-like.
        Time array in which to evaluate the dynamics.
    beta_s, beta_a, beta_aq : float.
        Infection rates for each of the three infected populations.
    gamma_s, gamma_a : float.
        Rate in which exposed individuals become infective with symptoms or
        asymptomatic, respectively.
    theta : float.
        Rate in which an asymptomatic individual goes into quarantine.
    mu : float.
        Rate in which infected individuals are removed from this group.
    Returns
    -------
    system_rhs : array-like.
        Returns evaluation of the right-hand side of the ODE system
    """
    # Unpack different population groups
    S, E, I_s, I_a, I_aq, R = x

    # Compute time derivatives
    dS_dt = -beta_s * S * I_s - beta_a * S * I_a - beta_aq * S * I_aq
    dE_dt = (
        beta_s * S * I_s
        + beta_a * S * I_a
        + beta_aq * S * I_aq
        - gamma_s * E
        - gamma_a * E
    )
    dI_s_dt = gamma_s * E - mu * I_s
    dI_a_dt = gamma_a * E - theta * I_a - mu * I_a
    dI_aq_dt = theta * I_a - mu * I_aq
    dR_dt = mu * (I_s + I_a + I_aq)

    # Return results as numpy array
    return np.array([dS_dt, dE_dt, dI_s_dt, dI_a_dt, dI_aq_dt, dR_dt])

To test that we can get the system up and running let's choose some simple initial condition and random set of parameter values. We'll then numerically integrate the system using scipy.integrate.odeint().

In [4]:
# Define time (in a.u.) for integration
t = np.linspace(0, 15, 500)

# Set initial conditions for population
init_infection = 1e-2
x0 = np.array(
    [
        1 - init_infection,  # S
        0,  # E
        init_infection / 3,  # I_s
        init_infection / 3,  # I_a
        init_infection / 3,  # I_aq
        0,  # R
    ]
)

# Arbitrary set of parameters
beta_s = 1
beta_a = 5
beta_aq = 1
gamma_s = 1
gamma_a = 3
theta = 1
mu = 0.5

# Pack parameters into tuple
args = (
    beta_s,
    beta_a,
    beta_aq,
    gamma_s,
    gamma_a,
    theta,
    mu,
)

# Numerically integrate system
x = sp.integrate.odeint(system_rhs, x0, t, args=args)

Let's now show the solution with an interactive Bokeh plot.

In [5]:
# Set colors for curves
colors = bokeh.palettes.d3["Category10"][5]

# Split populations
S, E, I_s, I_a, I_aq, R = x.transpose()

# Setup plot
p = bokeh.plotting.figure(
    plot_width=500,
    plot_height=300,
    x_axis_label="time (a.u.)",
    y_axis_label="fraction of population",
)

# Populate plot with glyphs
p.line(t, S, line_width=2, color=colors[0], legend_label="S")
p.line(t, E, line_width=2, color=colors[1], legend_label="E")
p.line(t, I_s, line_width=2, color=colors[2], legend_label="I_s")
p.line(t, I_a, line_width=2, color=colors[3], legend_label="I_a")
p.line(t, I_aq, line_width=2, color=colors[4], legend_label="I_aq")
p.line(t, R, line_width=2, color=colors[4], legend_label="R")

# Place legend
p.legend.location = "center_right"

# Show plot
bokeh.io.show(p)

The results follow the general SIR-like dynamics in which infected people initially increase exponentially and then decrease as the number of removed individual increases.

With a large parameter space let's turn our attention into making a set of interactive plots to explore this space.

Controlling parameters with panel

Now that we know that the system is behaving as we would intuitively expect, let's set the parameters such that we can control them via panel. This is a new library that allows easy generation of interactive widgets.

To have a better handling of the model in terms of intuitive parameters we will split some of the parameters into more sub-parameters that we can freely explore. For example the rate $\theta$ at which asymptomatic infected individuals go into quarantine will be written as a function of three parameters: The number of tests perform per unit time, the effectiveness of the test, and the probability of the person that tests positive complying with the quarantine order.

Let's first define all the widgets that we will need.

In [7]:
# Average infective contacts per day for symptomatic
contact_sym = pn.widgets.FloatSlider(
    name="infective contact symptomatics (βs)",
    start=0,
    end=2,
    step=0.05,
    value=0.1,
)
# Average infective contacts per day for asymptomatic
contact_asym = pn.widgets.FloatSlider(
    name="infective contacts asymptomatics (βa)",
    start=0,
    end=2,
    step=0.05,
    value=0.2,
)
# Average infective contacts per day for asymptomatic in quarantine
contact_asym_q = pn.widgets.FloatSlider(
    name="infective contacts quarantined (βaq)",
    start=0,
    end=2,
    step=0.05,
    value=0,
)
# Incubation period
incu_time = pn.widgets.FloatSlider(
    name="incubation period (Ɣ⁻¹) [days]", start=0, end=30, step=1, value=5
)
# Fraction of symptomatic patients
p_sym = pn.widgets.FloatSlider(
    name="symptomatic fraction", start=0, end=1, step=0.01, value=0.14
)
# Rate of testing (fraction of population per day)
k_test = pn.widgets.FloatSlider(
    name="rate of testing (fraction of population per day)",
    start=0,
    end=1,
    step=0.00001,
    value=0.0001,
)
# test efectiveness
delta_eff = pn.widgets.FloatSlider(
    name="test effectiveness (δeff)", start=0, end=1, step=0.05, value=0.70,
)
# test efectiveness
p_comp = pn.widgets.FloatSlider(
    name="probability of quarentine compliance",
    start=0,
    end=1,
    step=0.05,
    value=0.80,
)
# Time to recovery
recovery_time = pn.widgets.FloatSlider(
    name="time to recovery (µ⁻¹) [days]", start=7, end=30, step=1, value=15,
)
# Population size
log_pop = pn.widgets.FloatSlider(
    name="log pop. size", start=1, end=10, step=0.1, value=np.log10(3E8)
)
# Integration time
int_time = pn.widgets.IntSlider(
    name="integration time [days]", start=100, end=1000, step=1, value=100
)
# Plot cumulative infected people
infect_toggle = pn.widgets.Toggle(
    name='total infected', button_type='success'
)

Having defined these widgets we will now define a function that uses the values of these widgets, computes the dynamics, and returns an interactive plot of each of the populations.

In [8]:
@pn.depends(
    contact_sym.param.value,
    contact_asym.param.value,
    contact_asym_q.param.value,
    incu_time.param.value,
    p_sym.param.value,
    k_test.param.value,
    delta_eff.param.value,
    p_comp.param.value,
    recovery_time.param.value,
    log_pop.param.value,
    int_time.param.value,
    infect_toggle.param.value,
)
def SEIIIR_model(
    contact_sym,
    contact_asym,
    contact_asym_q,
    incu_time,
    p_sym,
    k_test,
    delta_eff,
    p_comp,
    recovery_time,
    log_pop,
    int_time,
    infect_toggle,
):
    """
    test
    """
    # Define time (in days) for integration
    t = np.linspace(0, int_time, 500)
    
    # Set parameters
    pop_size = 10**log_pop
    beta_s = contact_sym / pop_size
    beta_a = contact_asym / pop_size
    beta_aq = contact_asym_q / pop_size
    gamma = 1 / incu_time
    gamma_s = gamma * p_sym
    gamma_a = gamma * (1 - p_sym)
    theta = k_test * delta_eff * p_comp
    mu = 1 / recovery_time

    # Pack parameters into tuple
    args = (
        beta_s,
        beta_a,
        beta_aq,
        gamma_s,
        gamma_a,
        theta,
        mu,
    )
    
    # Set initial conditions for population
    init_infection = 1e-2
    x0 = np.array(
        [
            1 - init_infection,  # S
            0,  # E
            init_infection / 3,  # I_s
            init_infection / 3,  # I_a
            init_infection / 3,  # I_aq
            0,  # R
        ]
    ) * pop_size

    # Numerically integrate system
    x = sp.integrate.odeint(system_rhs, x0, t, args=args) / pop_size
    
    # Bokeh plot
    # Set colors for curves
    colors = bokeh.palettes.d3["Category10"][10]

    # Split populations
    S, E, I_s, I_a, I_aq, R = x.transpose()

    # Setup plot
    p = bokeh.plotting.figure(
        plot_width=500,
        plot_height=300,
        x_axis_label="time (days)",
        y_axis_label="fraction of population",
    )

    # Populate plot with glyphs
    p.line(t, S, line_width=2, color=colors[0], legend_label="S")
    p.line(t, E, line_width=2, color=colors[1], legend_label="E")
    if infect_toggle:
        p.line(
            t, I_s + I_a + I_aq,
            line_width=2,
            color=colors[6],
            legend_label="I"
        )
    else:
        p.line(t, I_s, line_width=2, color=colors[2], legend_label="I_s")
        p.line(t, I_a, line_width=2, color=colors[3], legend_label="I_a")
        p.line(t, I_aq, line_width=2, color=colors[4], legend_label="I_aq")
        
    p.line(t, R, line_width=2, color=colors[5], legend_label="R")

    # Place legend
    p.legend.location = "center_right"
    # Apply PBoC plot format
    covid19.viz.pboc_single(p)
    
    return(p)

Let's now arrange all parameter widgets into a dashboard to start exploring the parameters space.

In [9]:
pn.Column(
    pn.WidgetBox(
        pn.Row(
            pn.Column(contact_sym, contact_asym, contact_asym_q),
            pn.Column(incu_time, p_sym, recovery_time),
            pn.Column(k_test, delta_eff, p_comp)
        )
    ),
    pn.Row(
    SEIIIR_model, pn.WidgetBox(pn.Column(log_pop, int_time, infect_toggle))
    )
)
Out[9]:

Systematic exploration of the model

Maximum number of infected people.

An interesting quantity to explore under different scenarios what would be the maximum number of infected people at any point. In particular this is an interesting quantity to explore as a function of the number of tests available for the population. Let's again set an interactive plot, recycling some of the widgets generated for the previous plot.

In [10]:
# Define range of testing rates (fraction of pop. per day)
k_test_range = np.linspace(0, 1, 100)

# Establish function dependencies
@pn.depends(
    contact_sym.param.value,
    contact_asym.param.value,
    contact_asym_q.param.value,
    incu_time.param.value,
    p_sym.param.value,
    delta_eff.param.value,
    p_comp.param.value,
    recovery_time.param.value,
    log_pop.param.value,
)
def max_infection(
    contact_sym,
    contact_asym,
    contact_asym_q,
    incu_time,
    p_sym,
    delta_eff,
    p_comp,
    recovery_time,
    log_pop,
):
    # Define time (in days) for integration
    t = np.linspace(0, 500, 1000)

    # Set parameters
    pop_size = 10 ** log_pop
    beta_s = contact_sym / pop_size
    beta_a = contact_asym / pop_size
    beta_aq = contact_asym_q / pop_size
    gamma = 1 / incu_time
    gamma_s = gamma * p_sym
    gamma_a = gamma * (1 - p_sym)
    mu = 1 / recovery_time

    # Initialize array to save maximum number of infected cases
    I_max = np.zeros_like(k_test_range)

    # Loop through k_test values
    for i, k in enumerate(k_test_range):
        # Compute theta
        theta = k * delta_eff * p_comp
        # Pack parameters into tuple
        args = (
            beta_s,
            beta_a,
            beta_aq,
            gamma_s,
            gamma_a,
            theta,
            mu,
        )

        # Set initial conditions for population
        init_infection = 1e-2
        x0 = (
            np.array(
                [
                    1 - init_infection,  # S
                    0,  # E
                    init_infection / 3,  # I_s
                    init_infection / 3,  # I_a
                    init_infection / 3,  # I_aq
                    0,  # R
                ]
            )
            * pop_size
        )

        # Numerically integrate system
        x = sp.integrate.odeint(system_rhs, x0, t, args=args) / pop_size
        # Split populations
        S, E, I_s, I_a, I_aq, R = x.transpose()
        # Add number of infected and save it
        I_max[i] = np.max(I_s + I_a + I_aq)

    # Set Bokeh plot
    p = bokeh.plotting.figure(
        plot_width=500,
        plot_height=300,
        x_axis_label="fraction of population tested per day",
        y_axis_label="max. number of infections",
    )

    # Populate plot with glyphs
    p.line(
        k_test_range, I_max, line_width=2,
    )
    # Apply PBoC plot format
    covid19.viz.pboc_single(p)
    
    return(p)

As before, let's arrange the widgets into a dashboard.

In [11]:
pn.Column(
    pn.WidgetBox(
        pn.Row(
            pn.Column(contact_sym, contact_asym, contact_asym_q),
            pn.Column(incu_time, p_sym, recovery_time),
            pn.Column(delta_eff, p_comp, log_pop)
        )
    ),
    max_infection
)
Out[11]:

Parameter references

Incubation period

https://annals.org/aim/fullarticle/2762808/incubation-period-coronavirus-disease-2019-covid-19-from-publicly-reported

  • There were 181 confirmed cases with identifiable exposure and symptom onset windows to estimate the incubation period of COVID-19. The median incubation period was estimated to be 5.1 days (95% CI, 4.5 to 5.8 days), and 97.5% of those who develop symptoms will do so within 11.5 days (CI, 8.2 to 15.6 days) of infection. These estimates imply that, under conservative assumptions, 101 out of every 10 000 cases (99th percentile, 482) will develop symptoms after 14 days of active monitoring or quarantine.

Recovery time

https://els-jbs-prod-cdn.literatumonline.com/pb/assets/raw/Lancet/infographics/coronavirus/Coronavirus_MedianTimeline_Infographic-1584612208650.jpg